Here, we examine how lens size may affect lens transmission.

Data

First, we import the tidied data compiled from all studies that includes ASW names and ecological details.

#import species data for lens transmission
lens_data <- data.frame(read.csv("../Data/tidy data/lenses_compiled.csv", header=TRUE, na.strings=c("", "NA")))

Then we subset data to adults measured in this study only for analyses.

#adults-this study
lens_adults.new <- lens_data %>% 
  filter(stage %in% c("adult", "unknown")) %>%
  filter(source == "this_study")

Finally, we import the pruned tree.

#Import pruned tree
lens.tree <- read.tree(file = "../Data/tidy data/pruned_tree.txt")

Lens transmission vs. lens size

Here, we run PGLS regressions of t50 vs. log10 lens diameter and of %UVA vs. log10 lens diameter to test for relationships between transmission and path length.

Dataset: all adult species means

Here, we included only adult data from our study (we excluded published adult data because lens sizes from specimens used for transmission measurements were not reported). We include all species, whether or not lens transmission spectra appear to indicate the presence of a pigment. This was to test the null hypothesis that variation in transmission is due to different spectral properties in lenses across species. If we find a significant, highly correlated relationship between transmission and pathlength, it could indicate that all lenses have the same spectral properties, and differences in transmission may be driven only by differences in path length across species with different eye sizes.

t50 vs. path length

#PGLS models for t50 or %UVA vs. lens diameter among species

#### prep data ####

#Make list of taxa to drop (in tree but not in dataset)
drops <- setdiff(lens.tree$tip.label, as.character(lens_adults.new$ASW_names))

#Drop unwanted tips from phylogeny
adults.new.tree <- drop.tip(phy = lens.tree, tip = drops) 

#Make row names of the species the phylogeny tip labels
rownames(lens_adults.new) <- lens_adults.new$ASW_names

#check that names match in dataframe and tree
name.check(phy = adults.new.tree, data = lens_adults.new, data.names = lens_adults.new$ASW_names)

#use caper function to combine phylogeny and data into one object (this function also matches species names in tree and dataset)
lenssize.comp <- comparative.data(phy = adults.new.tree, data = lens_adults.new, 
                              names.col = ASW_names, vcv = TRUE, 
                              na.omit = FALSE, warn.dropped = TRUE)

#check for dropped tips or dropped species
lenssize.comp$dropped$tips #phylogeny
lenssize.comp$dropped$unmatched.rows #dataset

#### PGLS model for t50 vs. lens size ####

#model
pgls_size.t50 <- pgls(t50 ~ log10(approx_size), 
                  data = lenssize.comp, 
                  lambda = "ML", #uses Maximum Liklihood estimate of lambda
                  param.CI = 0.95)

After fitting the model, we can examine the diagnostic plots to check model assumptions, then look at the parameter estimates and the likelihood profile of lambda.

###check model assumptions ###

#diagnostic plots
par(mar = c(4,4,2,2))
par(mfrow = c(2, 2))
plot(pgls_size.t50, main = "PGLS model for t50 vs. log lens size")

par(mfrow = c(1, 1))

### model outputs ###

#print model output 
summary(pgls_size.t50)
## 
## Call:
## pgls(formula = t50 ~ log10(approx_size), data = lenssize.comp, 
##     lambda = "ML", param.CI = 0.95)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.9917 -0.9738  0.3471  1.7731  4.3893 
## 
## Branch length transformations:
## 
## kappa  [Fix]  : 1.000
## lambda [ ML]  : 0.926
##    lower bound : 0.000, p = 2.0995e-07
##    upper bound : 1.000, p = 0.0041166
##    95.0% CI   : (0.683, 0.991)
## delta  [Fix]  : 1.000
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)         345.530     15.171 22.7757   <2e-16 ***
## log10(approx_size)   13.355     11.383  1.1732   0.2441    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.942 on 83 degrees of freedom
## Multiple R-squared: 0.01631, Adjusted R-squared: 0.004462 
## F-statistic: 1.377 on 1 and 83 DF,  p-value: 0.2441
#Likelihood plot for Pagel's lambda from the PGLS model of eye diameter vs. the cuberoot of mass. Solid red line indicates estimate for lambda and broken red lines indcaite the 95% confidence interval
lambda.size.t50 <- pgls.profile(pgls_size.t50, "lambda")
plot(lambda.size.t50)

Finally, we can plot our data and the model fit. Note that the linear relationship between t50 and log lens diameter was not significant, indicating that observed variation in transmission among species is not driven solely by differences in path length.

#plot t50 vs. log lens size with fit
plot_t50size <- ggplot(lens_adults.new, aes(x = approx_size, y = t50, text = genus_species)) + 
  geom_point(size = 2, alpha = 0.8) + 
  scale_x_log10("Lens diameter (mm)") + 
  ylab("t50") + 
  theme_bw() +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  geom_abline(slope = coef(pgls_size.t50)[[2]], intercept = coef(pgls_size.t50)[[1]], linetype = "dashed")

#interactive plot
ggplotly(plot_t50size)

%UVA vs. path length

#### PGLS model for %UVA vs. lens size ####

#model
pgls_size.puva <- pgls(pUVA ~ log10(approx_size), 
                  data = lenssize.comp, 
                  lambda = "ML", #uses Maximum Liklihood estimate of lambda
                  param.CI = 0.95)

After fitting the model, we can examine the diagnostic plots to check model assumptions, then look at the parameter estimates and the likelihood profile of lambda.

###check model assumptions ###

#diagnostic plots
par(mar = c(4,4,2,2))
par(mfrow = c(2, 2))
plot(pgls_size.puva, main = "PGLS model for pUVA vs. log lens size")

par(mfrow = c(1, 1))

### model outputs ###

#print model output 
summary(pgls_size.puva)
## 
## Call:
## pgls(formula = pUVA ~ log10(approx_size), data = lenssize.comp, 
##     lambda = "ML", param.CI = 0.95)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.1903 -1.2283 -0.1703  0.5680  3.5408 
## 
## Branch length transformations:
## 
## kappa  [Fix]  : 1.000
## lambda [ ML]  : 0.839
##    lower bound : 0.000, p = 2.2887e-05
##    upper bound : 1.000, p = 0.00045649
##    95.0% CI   : (0.506, 0.973)
## delta  [Fix]  : 1.000
## 
## Coefficients:
##                    Estimate Std. Error t value  Pr(>|t|)    
## (Intercept)         55.6462    10.9241  5.0939 2.162e-06 ***
## log10(approx_size) -11.4285     9.1549 -1.2484    0.2154    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.446 on 83 degrees of freedom
## Multiple R-squared: 0.01843, Adjusted R-squared: 0.006604 
## F-statistic: 1.558 on 1 and 83 DF,  p-value: 0.2154
#Likelihood plot for Pagel's lambda from the PGLS model of eye diameter vs. the cuberoot of mass. Solid red line indicates estimate for lambda and broken red lines indcaite the 95% confidence interval
lambda.size.uva <- pgls.profile(pgls_size.puva, "lambda")
plot(lambda.size.uva)

Finally, we can plot our data and the model fit. Note that the linear relationship between %UVA and log lens diameter was not significant, indicating that observed variation in transmission among species is not driven solely by differences in path length.

#plot %UVA vs. log lens size with fit
plot_uvasize <- ggplot(lens_adults.new, aes(x = approx_size, y = pUVA, text = genus_species)) + 
  geom_point(size = 2, alpha = 0.8) + 
  scale_x_log10("Lens diameter (mm)") + 
  ylab("%UVA") + 
  theme_bw() +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  geom_abline(slope = coef(pgls_size.puva)[[2]], intercept = coef(pgls_size.puva)[[1]], linetype = "dashed")

#interactive plot
ggplotly(plot_uvasize)

Dataset: adult species means with unpigmented lenses

We only expect path length to show a strong relationship with UV transmission in unpigmented lenses, so we may want to look separately at just unpigmented ones. In our results, we considered a t50 ≤ 350nm to be an unpigmented lens (n = 27).

Our lenses ranged in t50 from 317 to 423. We can look at that plotted onto a histogram of all our t50 values to see where that falls in our data.

#Histogram of T50 values in our dataset
  plot_t50hist <-ggplot(lens_adults.new, aes(x = t50)) + 
  geom_histogram(color = "black", fill = "gray", binwidth = 5) +
  theme_bw() +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  xlab("t50") +
  geom_vline(aes(xintercept=350), col = "red", linetype = "dashed")
  

plot_t50hist

t50 vs. path length (unpigmented)

#add pigmentation category to dataframe
lens_adults.new <- lens_adults.new %>%
  mutate(pigment = ifelse(t50 <= 350, "unpigmented", "pigmented"))

#subset unpigmented species
lens_unpig <- lens_adults.new %>%
  filter(pigment=="unpigmented")

#PGLS models for t50 or %UVA vs. lens diameter among species

#### prep data ####

#Make list of taxa to drop (in tree but not in dataset)
drops <- setdiff(lens.tree$tip.label, as.character(lens_unpig$ASW_names))

#Drop unwanted tips from phylogeny
unpig.tree <- drop.tip(phy = lens.tree, tip = drops) 

#Make row names of the species the phylogeny tip labels
rownames(lens_unpig) <- lens_unpig$ASW_names

#check that names match in dataframe and tree
name.check(phy = unpig.tree, data = lens_unpig, data.names = lens_unpig$ASW_names)

#use caper function to combine phylogeny and data into one object (this function also matches species names in tree and dataset)
unpig.comp <- comparative.data(phy = unpig.tree, data = lens_unpig, 
                              names.col = ASW_names, vcv = TRUE, 
                              na.omit = FALSE, warn.dropped = TRUE)

#check for dropped tips or dropped species
unpig.comp$dropped$tips #phylogeny
unpig.comp$dropped$unmatched.rows #dataset

#### PGLS model for t50 vs. lens size ####

#model
pgls_unpig.t50 <- pgls(t50 ~ log10(approx_size), 
                  data = unpig.comp, 
                  lambda = "ML", #uses Maximum Liklihood estimate of lambda
                  param.CI = 0.95)

After fitting the model, we can examine the diagnostic plots to check model assumptions, then look at the parameter estimates and the likelihood profile of lambda.

###check model assumptions ###

#diagnostic plots
par(mar = c(4,4,2,2))
par(mfrow = c(2, 2))
plot(pgls_unpig.t50, main = "PGLS model for t50 vs. log lens size in unpigmented lenses")

par(mfrow = c(1, 1))

### model outputs ###

#print model output 
summary(pgls_unpig.t50)
## 
## Call:
## pgls(formula = t50 ~ log10(approx_size), data = unpig.comp, lambda = "ML", 
##     param.CI = 0.95)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.56134 -0.42041  0.05374  0.31628  1.46275 
## 
## Branch length transformations:
## 
## kappa  [Fix]  : 1.000
## lambda [ ML]  : 0.000
##    lower bound : 0.000, p = 1    
##    upper bound : 1.000, p = 0.028999
##    95.0% CI   : (NA, 0.962)
## delta  [Fix]  : 1.000
## 
## Coefficients:
##                    Estimate Std. Error  t value Pr(>|t|)    
## (Intercept)        333.7931     3.1406 106.2849   <2e-16 ***
## log10(approx_size)  10.1862     8.6350   1.1796   0.2488    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.6596 on 26 degrees of freedom
## Multiple R-squared: 0.0508,  Adjusted R-squared: 0.01429 
## F-statistic: 1.392 on 1 and 26 DF,  p-value: 0.2488
#Likelihood plot for Pagel's lambda from the PGLS model of eye diameter vs. the cuberoot of mass. Solid red line indicates estimate for lambda and broken red lines indcaite the 95% confidence interval
lambda.unpig.t50 <- pgls.profile(pgls_unpig.t50, "lambda")
plot(lambda.unpig.t50)

Finally, we can plot our data and the model fit. Note that the linear relationship between t50 and log lens diameter was not significant, indicating that lens transmission among species with unpigmented lenses is not correlated with path length.

#plot t50 vs. log lens size with fit
plot_t50unpig <- ggplot(lens_unpig, aes(x = approx_size, y = t50, text = genus_species)) + 
  geom_point(size = 2, alpha = 0.8) + 
  scale_x_log10("Lens diameter (mm)") + 
  ylab("t50") + 
  theme_bw() +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  geom_abline(slope = coef(pgls_unpig.t50)[[2]], intercept = coef(pgls_unpig.t50)[[1]], linetype = "dashed")

#interactive plot
ggplotly(plot_t50unpig)

%UVA vs. path length (unpigmented)

#### PGLS model for %UVA vs. lens size ####

#model
pgls_unpig.puva <- pgls(pUVA ~ log10(approx_size), 
                  data = unpig.comp, 
                  lambda = "ML", #uses Maximum Liklihood estimate of lambda
                  param.CI = 0.95)

After fitting the model, we can examine the diagnostic plots to check model assumptions, then look at the parameter estimates and the likelihood profile of lambda.

###check model assumptions ###

#diagnostic plots
par(mar = c(4,4,2,2))
par(mfrow = c(2, 2))
plot(pgls_unpig.puva, main = "PGLS model for pUVA vs. log lens size in unpigmented lenses")

par(mfrow = c(1, 1))

### model outputs ###

#print model output 
summary(pgls_unpig.puva)
## 
## Call:
## pgls(formula = pUVA ~ log10(approx_size), data = unpig.comp, 
##     lambda = "ML", param.CI = 0.95)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.3417 -0.2151  0.0785  0.1825  1.2231 
## 
## Branch length transformations:
## 
## kappa  [Fix]  : 1.000
## lambda [ ML]  : 0.000
##    lower bound : 0.000, p = 1    
##    upper bound : 1.000, p = 0.011065
##    95.0% CI   : (NA, 0.863)
## delta  [Fix]  : 1.000
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        63.35822    2.66112 23.8088   <2e-16 ***
## log10(approx_size) -0.48154    7.31682 -0.0658    0.948    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.5589 on 26 degrees of freedom
## Multiple R-squared: 0.0001666,   Adjusted R-squared: -0.03829 
## F-statistic: 0.004331 on 1 and 26 DF,  p-value: 0.948
#Likelihood plot for Pagel's lambda from the PGLS model of eye diameter vs. the cuberoot of mass. Solid red line indicates estimate for lambda and broken red lines indcaite the 95% confidence interval
lambda.unpig.uva <- pgls.profile(pgls_unpig.puva, "lambda")
plot(lambda.unpig.uva)

Finally, we can plot our data and the model fit. Note that the linear relationship between t50 and log lens diameter was not significant, indicating that lens transmission among species with unpigmented lenses is not correlated with path length.

#plot %UVA vs. log lens size with fit
plot_uva.unpig <- ggplot(lens_unpig, aes(x = approx_size, y = pUVA, text = genus_species)) + 
  geom_point(size = 2, alpha = 0.8) + 
  scale_x_log10("Lens diameter (mm)") + 
  ylab("%UVA") + 
  theme_bw() +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  geom_abline(slope = coef(pgls_unpig.puva)[[2]], intercept = coef(pgls_unpig.puva)[[1]], linetype = "dashed")

#interactive plot
ggplotly(plot_uva.unpig)

The lack of relationship between transmission and lens size across species is good evidence that differences we observe across species are due to differences in the spectral properties of the lens, even among species with highly transparent lenses. It also justifies why we do not account for lens diameter in our models of lens transmission vs. ecology (which is impossible, as the published data we include did not have measurements for lens diameter).

Ecology vs. lens size

Here we look at whether lens size varies across ecological traits in our dataset.

# create vectors of colors for habitat
col_hab <- c("scansorial" = "#009E73",
              "other"  = "#636363")

#create vector of colors for activity period
col_act <- c("nondiurnal" = "#636363",
             "diurnal" = "#FFA010")

#boxplot of lens size across scansorial/not
plot_scans <- ggplot(lens_adults.new, 
                aes(x = hab, y = approx_size, fill = hab)) +
  geom_violin(trim=FALSE, alpha = 0.6)+
  scale_fill_manual(values = col_hab, name = " ") +
  scale_x_discrete(breaks = c("scansorial","other"),
                   limits = c("scansorial","other"),
                   labels = c("scansorial","non-scansorial")) +
  labs(title=NULL,x=" ", y = "Lens diameter (mm)")+
  geom_boxplot(width=0.15, fill="white", outlier.alpha = 0)+
  geom_jitter(shape = 19, size = 1, alpha = 0.7, position = position_jitter(0.1)) + 
  theme_classic() +
  theme(legend.position = "none")

#boxplot of lens size across diurnal/not
plot_diur <- ggplot(lens_adults.new, 
                aes(x = act, y = approx_size, fill = act)) +
  geom_violin(trim=FALSE, alpha = 0.6)+
  scale_fill_manual(values = col_act, name = " ") +
  scale_x_discrete(breaks = c("diurnal","nondiurnal"),
                   limits = c("diurnal","nondiurnal"),
                   labels = c("diurnal","non-diurnal")) +
  labs(title=NULL,x=" ", y = "Lens diameter (mm)")+
  geom_boxplot(width=0.15, fill="white", outlier.alpha = 0)+
  geom_jitter(shape = 19, size = 1, alpha = 0.7, position = position_jitter(0.1)) + 
  theme_classic() +
  theme(legend.position = "none")

#make figure
plot_grid(plot_scans, plot_diur + ylab(""),
          align = 'vh', 
          labels = c("A", "B"),
          hjust = -0.5, 
          vjust = 1.8, 
          nrow = 1)

#export fig
pdf("../Figures/FigS7.pdf", width = 6, height = 4)
plot_grid(plot_scans, plot_diur + ylab(""),
          align = 'vh', 
          labels = c("A", "B"),
          hjust = -0.5, 
          vjust = 1.8, 
          nrow = 1)
dev.off()
## quartz_off_screen 
##                 2

We can also test whether species with different ecological traits differ in lens size.

#PGLS model for lens size vs. scansoriality -----

#model
pgls_hab <- pgls(approx_size ~ hab, 
                  data = lenssize.comp, 
                  lambda = "ML", #uses Maximum Liklihood estimate of lambda
                  param.CI = 0.95)

#diagnostic plots
par(mar = c(4,4,2,2))
par(mfrow = c(2, 2))
plot(pgls_hab, main = "PGLS model for size vs scansoriality")

par(mfrow = c(1, 1))

#print model output 
summary(pgls_hab)
## 
## Call:
## pgls(formula = approx_size ~ hab, data = lenssize.comp, lambda = "ML", 
##     param.CI = 0.95)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.23999 -0.03812  0.02549  0.10064  0.25369 
## 
## Branch length transformations:
## 
## kappa  [Fix]  : 1.000
## lambda [ ML]  : 0.769
##    lower bound : 0.000, p = 0.0076294
##    upper bound : 1.000, p = 0.0075615
##    95.0% CI   : (0.220, 0.974)
## delta  [Fix]  : 1.000
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)  
## (Intercept)    1.94056    0.75853  2.5583  0.01233 *
## habscansorial  0.33297    0.41720  0.7981  0.42709  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1058 on 83 degrees of freedom
## Multiple R-squared: 0.007616,    Adjusted R-squared: -0.004341 
## F-statistic: 0.637 on 1 and 83 DF,  p-value: 0.4271
#PGLS model for lens size vs. diurnality -----

#model
pgls_act <- pgls(approx_size ~ act, 
                  data = lenssize.comp, 
                  lambda = "ML", #uses Maximum Liklihood estimate of lambda
                  param.CI = 0.95)

#diagnostic plots
par(mar = c(4,4,2,2))
par(mfrow = c(2, 2))
plot(pgls_act, main = "PGLS model for size vs diurnality")

par(mfrow = c(1, 1))

#print model output 
summary(pgls_act)
## 
## Call:
## pgls(formula = approx_size ~ act, data = lenssize.comp, lambda = "ML", 
##     param.CI = 0.95)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.24674 -0.02451  0.01147  0.06700  0.20775 
## 
## Branch length transformations:
## 
## kappa  [Fix]  : 1.000
## lambda [ ML]  : 0.641
##    lower bound : 0.000, p = 0.0077929
##    upper bound : 1.000, p = 0.0010038
##    95.0% CI   : (0.120, 0.939)
## delta  [Fix]  : 1.000
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)  
## (Intercept)    0.93010    0.79779  1.1658  0.24702  
## actnondiurnal  1.09399    0.49781  2.1976  0.03076 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.09542 on 83 degrees of freedom
## Multiple R-squared: 0.05499, Adjusted R-squared: 0.0436 
## F-statistic:  4.83 on 1 and 83 DF,  p-value: 0.03076

These models indicate that lens size does not differ significantly among scansorial vs. nonscansorial species, but that nondiurnal species have significantly larger lenses than diurnal species in this dataset.

Supplemental Figure

#make panels to figures
fig.a <- plot_t50size +
  xlab("") +
  ylab("t50 (nm)") + 
  theme_classic()

fig.b <- plot_uvasize +
  xlab("")+ 
  theme_classic()

fig.c <- plot_t50unpig +
  ylab("t50 (nm)")+ 
  theme_classic()

fig.d <- plot_uva.unpig+ 
  theme_classic()

fig.e <- plot_scans

fig.f <- plot_diur

# arrange panels in figure
plots <- plot_grid(fig.a, fig.b, fig.c, fig.d, fig.e, fig.f,
           align = 'vh', 
           labels = c("A", "B", "C", "D", "E", "F"), #panel labels for figure
           hjust = -.5,
           vjust = 1,
           nrow = 3) #number of rows in grids

#view figure
plot_grid(plots)

#export figure
pdf("../Figures/Fig-S6.pdf", width = 6, height = 7.5)
plot_grid(plots)
dev.off()
## quartz_off_screen 
##                 2